{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://cs229.stanford.edu/notes/cs229-notes3.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Understanding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functional margin ($\\hat \\gamma$) and geometric margin ($\\gamma$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definition for the functional margin of data point $x^{(i)}$:\n",
    "\n",
    "$$\\hat \\gamma^{(i)} = y^{(i)}(w^T x^{(i)} + b)$$\n",
    "\n",
    "Given $y^{(i)}$ is $\\pm 1$ for the positive and negative examples, $\\hat \\gamma^{(i)}$ is actually equivalent to\n",
    "\n",
    "$$| w^T x^{(i)} + b |$$\n",
    "\n",
    "which is a scaled distance of $x$ to the line (verified by geometry)\n",
    "\n",
    "$$w^T x^{(i)} + b = 0 $$ \n",
    "\n",
    "and the scaler is $\\left\\|w\\right\\|$. \n",
    "\n",
    "If we let $x_0$ be the vector that is perpendicular (hence $x_0$ is more specific to the $B$ shown in the lecture note) to the line, and also exactly on the line. Then $b$ can be written as \n",
    "\n",
    "$$b=-\\left\\|w\\right\\| \\left\\|x_0 \\right\\|$$\n",
    "\n",
    "$\\hat \\gamma$ can also be written in terms of $x_0$. To sum up:\n",
    "\n",
    "\\begin{align*}\n",
    "\\hat \\gamma^{(i)} &= y^{(i)}(w^T x^{(i)} + b)\\\\\n",
    "            &= |w^T x^{(i)} + b| \\\\\n",
    "            &= y^{(i)}w^T(x^{(i)} - x_0)\n",
    "\\end{align*}\n",
    "\n",
    "Note when $\\left\\|w\\right\\| = 1$. $\\hat \\gamma$ becomes the unscaled distance, i.e. the **geometric margin** ($\\gamma$, no hat).\n",
    "\n",
    "\n",
    "In the context of a **training set** with multiple data points ($x^{(1)}$, $x^{(2)}$, $\\cdots$, $x^{(m)}$), we define the functional margin of the training set ($\\hat \\Gamma$) to be \n",
    "\n",
    "$$\\hat \\Gamma = \\mathrm{min}_{i=i,\\cdots,m} \\hat \\gamma^{(i)}$$\n",
    "\n",
    "Note in the lecture notes, they use $\\hat \\gamma$ to represent the functional margin of a training set, but I feel it rather confusing as $\\hat \\gamma{(i)}$ has already been used for the functional margin of individual data point, so I choose to use upper case $\\hat \\Gamma$, which feels less confusing. Similarly, the geometric margin of a training set would be $\\Gamma$ without the hat sign."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SVM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After some reasoning: the SVM becomes an optimization problem for the following optimal margin classifier:\n",
    "\n",
    "$$\n",
    "\\mathrm{min}_{\\gamma, w, b} \\frac{1}{2}\\left\\|w\\right\\|^2\n",
    "$$\n",
    "Subject to:\n",
    "$$y^{(i)}(w^Tx^{(i)} + b) \\ge 1, i = 1, 2, \\cdots, m$$\n",
    "\n",
    "Here we enforced the functional margin to be exactly one ($\\hat \\gamma = 1$) after after scaling $w$ and $b$. Meanwhile, we want to minimize $\\left\\|w\\right\\|$.\n",
    "\n",
    "**Important to think through**: this is equivalent to maximizing the geometric margin when enforcing $\\left\\|w\\right\\| = 1$, given which, $y^{(i)}(w^Tx^{(i)} + b)$ can be interpreted as exactly the distance to the decision boundary.\n",
    "\n",
    "Recall \n",
    "\n",
    "$$\\hat \\gamma^{(i)} = y^{(i)}(w^T x^{(i)} + b)$$\n",
    "\n",
    "It's just about either\n",
    "\n",
    "1. fix $f(w) = \\left\\|w\\right\\| = 1$ involved in the right hand, maximize the left hand (the original idea behind maximum margin)\n",
    "1. fix left hand (because scaling $w$ and $b$ won't affect decision boundry line), minimize $f(w) = \\left\\|w\\right\\|$ which is involved in the right hand\n",
    "\n",
    "The two approaches have the same effect as to which decision boundry will be selected after the optimization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying Lagrange duality (non-trivial, esp. for the linearly unseparable cases), the optimization problem becomes\n",
    "\n",
    "$$\n",
    "\\mathrm{max}_{\\alpha} W(\\alpha) = \\sum_{i=1}^{m} \\alpha_i - \\frac{1}{2} \\sum_{i,j=1}^{m}y^{(i)}y^{(j)} \\alpha_i \\alpha_j x_i^T x_j\n",
    "$$\n",
    "Subject to:\n",
    "$$0 \\le \\alpha_i \\le C, i = 1, 2, \\cdots, m$$\n",
    "and\n",
    "$$\\sum_{i=1}^{m} \\alpha_i y^{(i)} = 0$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sequential Minimal optimization (SMO) algorithm is a relatively recent way to solve the above optimization problem. There are other ways, too. See John Platt's paper: https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/smo-book.pdf. SMO is also implemented in [libsvm](https://en.wikipedia.org/wiki/LIBSVM)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Follow http://cs229.stanford.edu/materials/smo.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "e.g. while optimizaing $\\alpha_1, \\alpha_2$, fix the other $\\alpha$s,\n",
    "\n",
    "\\begin{align*} \n",
    "W(\\alpha) &= \\alpha_1 + \\alpha_2 + \\sum_{i=3}^{m} \\alpha_i + \\frac{1}{2} [\\\\\n",
    "          &y_1^2 \\alpha_1^2 x_1^T x_1 + \\\\\n",
    "          &y_2^2 \\alpha_2^2 x_2^T x_2 + \\\\\n",
    "          &y_1 y_2 \\alpha_1 \\alpha_2 x_1^T x_2 + \\\\\n",
    "          &y_2 y_1 \\alpha_2 \\alpha_1 x_2^T x_1 + \\\\\n",
    "          &\\sum_{j=3}^{m} y_1 y_j \\alpha_1 \\alpha_j x_1^T x_j + \\\\\n",
    "          &\\sum_{j=3}^{m} y_2 y_j \\alpha_2 \\alpha_j x_2^T x_j + \\\\\n",
    "          &\\sum_{i=3}^{m} y_i y_1 \\alpha_i \\alpha_1 x_i^T x_1 + \\\\\n",
    "          &\\sum_{i=3}^{m} y_i y_2 \\alpha_i \\alpha_2 x_i^T x_2 + \\\\\n",
    "          &\\sum_{i=3,j=3}^{m} y_i y_j \\alpha_i \\alpha_j x_i^T x_j \\\\\n",
    "]\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# toy dataset\n",
    "Xs = np.array([\n",
    "        [0.5, 0.5],\n",
    "        [1, 1],\n",
    "        [1.5, 2],\n",
    "        [2, 1]\n",
    "    ])\n",
    "ys = [-1, -1, 1, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x10fdb7438>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEkhJREFUeJzt3W+MneV95vHvFWzQItOQ4qk3NZCBFYTa2dKSaUiXaNfd\naGtDUdkiVhuCgtZqZaFNqlTqC6JIBYm8qKoqVRTRxLKocaMgeJFYDUEUb6Rt4i4JUcYRwR5bQV5D\nwEDtAVYpfySD5d++OMdmGGZ8jmfOzBnf8/1I1sx5nlvnvjB+Lt++z3PmpKqQJLXlfcMOIEkaPMtd\nkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1KAVw5p49erVNTo6OqzpJemstGfPnper\naqTXuKGV++joKOPj48OaXpLOSkl+0c84t2UkqUGWuyQ1yHKXpAZZ7pLUIMtdkhrUs9yTXJLkn5Ls\nTzKR5PMzjEmSryY5mOSpJNcsTFxJfTn2Guz5e/je3Z2vx14bdiItsn5uhTwO/HlV/TTJBcCeJN+r\nqv1TxlwPXNH9dS3w9e5XSYvtFz+CB26BOgFvvwkrz4ddX4TbvgUf+t1hp9Mi6blyr6qXquqn3e9f\nAw4Aa6cNuwn4RnU8AVyY5IMDTyvp9I691in2t17vFDt0vr71euf4sdeHm0+L5oz23JOMAr8N/Hja\nqbXA81MeH+a9fwFIWmj7dnZW7DOpEzCxc3HzaGj6Lvckq4BvA39WVf86l8mSbEkynmR8cnJyLk8h\n6XRePfTOin26t9+EVw4tbh4NTV/lnmQlnWJ/oKpm+qv/BeCSKY8v7h57l6raVlVjVTU2MtLzRyNI\nOlO/enlnj30mK8+Hiy5f3Dwamn7ulgnwd8CBqvqbWYY9DNzevWvm48Avq+qlAeaU1I+P3AyZ5bLO\n+2D9zYubR0PTz90y1wGfAfYmebJ77IvApQBVtRV4FLgBOAi8CWwefFRJPZ13QeeumOl3y+R9nePn\nrRp2Qi2SnuVeVf8HSI8xBXx2UKEkzcOHfhf+/OedF09fOdTZill/s8W+zAztR/5KWkDnrYJrbh92\nCg2RP35AkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUu\nSQ2y3CWpQZa7JDXIcpekBlnuktSgfj5DdXuSo0n2zXL+/Um+m+RnSSaS+BF7kjRk/azcdwCbTnP+\ns8D+qroa2AB8Ocm5848mSZqrnuVeVbuBV083BLggSYBV3bHHBxNPkjQXg/gM1XuBh4EXgQuA/15V\nJwbwvJKkORrEC6obgSeBXwd+C7g3ya/MNDDJliTjScYnJycHMLUkaSaDKPfNwM7qOAg8A1w108Cq\n2lZVY1U1NjIyMoCpJUkzGUS5Pwd8EiDJGuDDwKEBPK8kaY567rkneZDOXTCrkxwG7gZWAlTVVuBL\nwI4ke4EAd1bVywuWWJLUU89yr6pbe5x/Efj9gSWSJM2b71CVpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ\n5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnu\nktSgnuWeZHuSo0n2nWbMhiRPJplI8oPBRpQknal+Vu47gE2znUxyIfA14A+raj3w3wYTTZI0Vz3L\nvap2A6+eZsingZ1V9Vx3/NEBZZMkzdEg9tyvBD6Q5PtJ9iS5fbaBSbYkGU8yPjk5OYCpJUkzGUS5\nrwA+CvwBsBH4iyRXzjSwqrZV1VhVjY2MjAxgaknSTFYM4DkOA69U1RvAG0l2A1cDTw/guSVJczCI\nlft3gE8kWZHkfOBa4MAAnleSNEc9V+5JHgQ2AKuTHAbuBlYCVNXWqjqQ5DHgKeAEcF9VzXrbpCRp\n4fUs96q6tY8xfw389UASSZLmzXeoSlKDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpk\nuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoN6lnuS7UmOJjntR+cl+Z0k\nx5PcMrh4kqS56GflvgPYdLoBSc4B/gr4XwPIJEmap57lXlW7gVd7DPtT4NvA0UGEkiTNz7z33JOs\nBf4I+HofY7ckGU8yPjk5Od+pJUmzGMQLql8B7qyqE70GVtW2qhqrqrGRkZEBTC1JmsmKATzHGPBQ\nEoDVwA1JjlfVPwzguSVJczDvcq+qy05+n2QH8IjFLknD1bPckzwIbABWJzkM3A2sBKiqrQuaTpI0\nJz3Lvapu7ffJqup/zCuNJGkgfIeqJDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGW\nuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNahnuSfZnuRokn2znL8tyVNJ\n9ib5YZKrBx9TknQm+lm57wA2neb8M8B/qqp/D3wJ2DaAXJKkeejnM1R3Jxk9zfkfTnn4BHDx/GNJ\nkuZj0Hvufwz842wnk2xJMp5kfHJycsBTS5JOGli5J/k9OuV+52xjqmpbVY1V1djIyMigppYkTdNz\nW6YfSX4TuA+4vqpeGcRzSpLmbt4r9ySXAjuBz1TV0/OPJEmar54r9yQPAhuA1UkOA3cDKwGqaitw\nF3AR8LUkAMeramyhAkuSeuvnbplbe5z/E+BPBpZIkjRvvkNVkhpkuUtSgyx3SWqQ5S5JDbLcJalB\nlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy13zUlU8\ntu9fqKq+jkvL1rHXYM/fw/fu7nw99tqCTtez3JNsT3I0yb5ZzifJV5McTPJUkmsGH1NL1a6JI9zx\nzT3c88j+U0VeVdzzyH7u+OYedk0cGXJCaQn4xY/gy1fBY1+Ax7/S+frlqzrHF0g/K/cdwKbTnL8e\nuKL7awvw9fnH0tli4/o1bL5ulPsff/ZUwd/zyH7uf/xZNl83ysb1a4YdURquY6/BA7fAW6/D2292\njr39ZufxA7fAsdcXZNp+PkN1d5LR0wy5CfhGdZZtTyS5MMkHq+qlAWXUEpaEu25cB8D9jz/L/Y8/\nC8Dm60a568Z1dD80XVq+9u2EOjHzuToBEzvhmtsHPu0g9tzXAs9PeXy4e+w9kmxJMp5kfHJycgBT\naymYWvAnWexS16uH3lmxT/f2m/DKoQWZdlFfUK2qbVU1VlVjIyMjizm1FtDJrZippu7BS8var14O\nK8+f+dzK8+Giyxdk2kGU+wvAJVMeX9w9pmVg+h77M395w3v24KVl7SM3Q2ap2rwP1t+8INMOotwf\nBm7v3jXzceCX7rcvH7smjpwq9pNbMXfduO5UwXu3jJa98y6A274F5656ZwW/8vzO49u+BeetWpBp\n02tlleRBYAOwGjgC3A2sBKiqrelsrN5L546aN4HNVTXea+KxsbEaH+85TEtcVbFr4ggb16951x77\nbMelZevY650XT1851NmKWX/znIo9yZ6qGus5blj/bLbcJenM9VvuvkNVkhpkuUtSgyx3SWqQ5S5J\nDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQg\ny12SGtRXuSfZlOTnSQ4m+cIM59+f5LtJfpZkIsnmwUeVJPWrZ7knOQf4W+B6YB1wa5J104Z9Fthf\nVVfT+bzVLyc5d8BZJUl96mfl/jHgYFUdqqq3gIeAm6aNKeCC7odlrwJeBY4PNKkkqW/9lPta4Pkp\njw93j011L/AbwIvAXuDzVXVi+hMl2ZJkPMn45OTkHCNLknoZ1AuqG4EngV8Hfgu4N8mvTB9UVduq\naqyqxkZGRgY0tSRpun7K/QXgkimPL+4em2ozsLM6DgLPAFcNJqIk6Uz1U+4/Aa5Icln3RdJPAQ9P\nG/Mc8EmAJGuADwOHBhlUktS/Fb0GVNXxJJ8DdgHnANuraiLJHd3zW4EvATuS7AUC3FlVLy9gbknS\nafQsd4CqehR4dNqxrVO+fxH4/cFGkyTNle9QlaQGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWp\nQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAb1Ve5JNiX5eZKD\nSb4wy5gNSZ5MMpHkB4ONKUk6Ez0/Zi/JOcDfAv8FOAz8JMnDVbV/ypgLga8Bm6rquSS/tlCBJUm9\n9bNy/xhwsKoOVdVbwEPATdPGfBrYWVXPAVTV0cHGlCSdiX7KfS3w/JTHh7vHproS+ECS7yfZk+T2\nmZ4oyZYk40nGJycn55ZYktTToF5QXQF8FPgDYCPwF0munD6oqrZV1VhVjY2MjAxoaknSdD333IEX\ngEumPL64e2yqw8ArVfUG8EaS3cDVwNMDSSlJOiP9rNx/AlyR5LIk5wKfAh6eNuY7wCeSrEhyPnAt\ncGCwUSVJ/eq5cq+q40k+B+wCzgG2V9VEkju657dW1YEkjwFPASeA+6pq30IGlyTNLlU1lInHxsZq\nfHx8KHNL0tkqyZ6qGus1zneoSlKDLHdJapDlLkkNOqvKvap4bN+/MP11gtmOS9JydVaV+66JI9zx\nzT3c88j+U0VeVdzzyH7u+OYedk0cGXJCSVoa+nkT05Kxcf0aNl83yv2PPwvAXTeu455H9nP/48+y\n+bpRNq5fM9yAkrREnFXlnoS7blwHwP2PP3uq5DdfN8pdN64jyRDTSdLScVZty8C7C/4ki12S3u2s\nK/eTe+xTTd2DlySdZeV+sthP7rE/85c3nNqDt+Al6R1n1Z77rokjp4r95FbM1D34ay+7iE0f+bdD\nTilJw3dW/WyZqmLXxBE2rl/zrj322Y5LUmv6/dkyZ9XKPcmMK/PZjkvScnVW7blLkvpjuUtSgyx3\nSWqQ5S5JDbLcJalBlrskNWho97knmQR+MZTJZ7caeHnYIU5jqeeDpZ/RfPOz1PPB0s8433wfqqqR\nXoOGVu5LUZLxft4cMCxLPR8s/Yzmm5+lng+WfsbFyue2jCQ1yHKXpAZZ7u+2bdgBeljq+WDpZzTf\n/Cz1fLD0My5KPvfcJalBrtwlqUHLstyTbEry8yQHk3xhljEbkjyZZCLJD5ZSviTvT/LdJD/r5tu8\nyPm2JzmaZN8s55Pkq938TyW5Zonlu62ba2+SHya5einlmzLud5IcT3LLYmWbMnfPjEO+Rnr9Px72\nNXJJkn9Ksr87/+dnGLOw10lVLatfwDnA/wUuB84FfgasmzbmQmA/cGn38a8tsXxfBP6q+/0I8Cpw\n7iJm/I/ANcC+Wc7fAPwjEODjwI8X+f9xr3z/AfhA9/vrl1q+KX8O/jfwKHDLYubr8/dwaNdIn/mG\nfY18ELim+/0FwNMzXMcLep0sx5X7x4CDVXWoqt4CHgJumjbm08DOqnoOoKqOLrF8BVyQzieTrKLz\nB/f4YgWsqt3dOWdzE/CN6ngCuDDJBxcnXe98VfXDqvp/3YdPABcvSrB35u/1+wfwp8C3gcX8s3dK\nHxmHeY30k2/Y18hLVfXT7vevAQeAtdOGLeh1shzLfS3w/JTHh3nvb/qVwAeSfD/JniS3L1q6/vLd\nC/wG8CKwF/h8VZ1YnHh96ee/Yan4YzqrpyUjyVrgj4CvDzvLaQzzGunHkrlGkowCvw38eNqpBb1O\nzqpPYlpEK4CPAp8E/g3woyRPVNXTw411ykbgSeA/A/8O+F6Sf66qfx1urLNLkt+jU+6fGHaWab4C\n3FlVJ5bwx0Z6jfQhySo6/wL7s8Weezmu3F8ALpny+OLusakOA7uq6o2qehnYDSzWi2795NtM55/E\nVVUHgWeAqxYpXz/6+W8YqiS/CdwH3FRVrww7zzRjwENJngVuAb6W5L8ON9J7DPMa6cfQr5EkK+kU\n+wNVtXOGIQt6nSzHcv8JcEWSy5KcC3wKeHjamO8An0iyIsn5wLV09syWSr7n6KyYSLIG+DBwaJHy\n9eNh4Pbu3QAfB35ZVS8NO9RJSS4FdgKfWUIrzVOq6rKqGq2qUeBbwP+sqn8YcqzphnmN9GOo10h3\nr//vgANV9TezDFvQ62TZbctU1fEknwN20bkjYXtVTSS5o3t+a1UdSPIY8BRwArivqk5729pi5gO+\nBOxIspfOK+13dldPiyLJg8AGYHWSw8DdwMop+R6lcyfAQeBNOquoRdNHvruAi+isiAGO1yL+oKk+\n8g1dr4zDvEb6yceQrxHgOuAzwN4kT3aPfRG4dErGBb1OfIeqJDVoOW7LSFLzLHdJapDlLkkNstwl\nqUGWuyQ1yHKXpAZZ7pLUIMtdkhr0/wEGgp+OJcRkFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c52b710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(*Xs[:2].T, marker='x', s=50)\n",
    "plt.scatter(*Xs[2:].T, marker='o', s=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calc_f(xk, alphas, Xs, ys, b):\n",
    "    return sum(alpha * yi * xi.dot(xk) for (alpha, xi, yi) in zip(alphas, Xs, ys)) + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calc_E(xk, yk, alphas, Xs, ys, b):\n",
    "    return calc_f(xk, alphas, Xs, ys, b) - yk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calc_L(yi, yj, ai, aj, C):\n",
    "    if yi != yj:\n",
    "        return max(0, aj - ai)\n",
    "    elif yi == yj:\n",
    "        return max(0, ai + aj - C)\n",
    "    else:\n",
    "        raise\n",
    "        \n",
    "def calc_H(ai, aj, yi, yj, C):\n",
    "    if yi != yj:\n",
    "        return min(C, aj - ai + C)\n",
    "    elif yi == yj:\n",
    "        return min(C, ai + aj)\n",
    "    else:\n",
    "        raise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def update_aj(aj, yj, Ei, Ej, eta, H, L):\n",
    "    \"\"\"Eq. 12 and 15\"\"\"\n",
    "    aj = aj - yj * (Ei - Ej) / eta\n",
    "    if aj > H:\n",
    "        return H\n",
    "    elif aj < L:\n",
    "        return L\n",
    "    else:\n",
    "        return aj\n",
    "\n",
    "def update_ai(ai, yi, yj, aj_old, aj):\n",
    "    \"\"\"Eq. 16\"\"\"\n",
    "    return ai + yi * yj * (aj_old - aj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calc_b(b1, b2, ai, aj, C):\n",
    "    if 0 < ai < C:\n",
    "        return b1\n",
    "    elif 0 < aj < C:\n",
    "        return b2\n",
    "    else:\n",
    "        return (b1 + b2) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# the toy example is linearly separable, so C could be arbitrarily large\n",
    "C = 1e8\n",
    "tol = 0.001\n",
    "max_passes = 10000\n",
    "\n",
    "m = Xs.shape[0]\n",
    "alphas = np.zeros(m)\n",
    "b = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "passes = 0\n",
    "while passes < max_passes:\n",
    "#     print(passes, end=',')\n",
    "    num_changed_alphas = 0\n",
    "    for i in range(m):\n",
    "        ai = alphas[i]\n",
    "        xi = Xs[i]\n",
    "        yi = ys[i]\n",
    "        Ei = calc_E(xi, yi, alphas, Xs, ys, b)\n",
    "        if (yi * Ei < -tol and ai < C) or (yi * Ei > tol and ai > 0):\n",
    "            j = np.random.choice([_ for _ in range(m) if _ != i])\n",
    "            aj = alphas[j]\n",
    "            xj = Xs[j]\n",
    "            yj = ys[j]\n",
    "            Ej = calc_E(xj, yj, alphas, Xs, ys, b)\n",
    "            \n",
    "            ai_old = ai\n",
    "            aj_old = aj\n",
    "            \n",
    "            L = calc_L(ai, aj, yi, yj, C)\n",
    "            H = calc_H(ai, aj, yi, yj, C)\n",
    "            if L == H:\n",
    "                continue\n",
    "                \n",
    "            eta = 2 * xi.dot(xj) - xi.dot(xi) - xj.dot(xj)\n",
    "            if eta >= 0:\n",
    "                continue\n",
    "                \n",
    "            aj = update_aj(aj, yj, Ei, Ej, eta, H, L)\n",
    "            alphas[j] = aj\n",
    "            if np.abs(aj - aj_old) < 1e-5:\n",
    "                continue\n",
    "            \n",
    "            ai = update_ai(ai, yi, yj, aj_old, aj)\n",
    "            alphas[i] = ai\n",
    "            \n",
    "            # Eq. 17\n",
    "            b1 = b - Ei - yi * (ai - ai_old) * xi.dot(xi) - yj * (aj - aj_old) * xi.dot(xj)\n",
    "            # Eq. 18\n",
    "            b2 = b - Ej - yi * (ai - ai_old) * xi.dot(xj) - yj * (aj - aj_old) * xj.dot(xj)\n",
    "            b = calc_b(b1, b2, ai, aj, C)\n",
    "            num_changed_alphas += 1\n",
    "    if num_changed_alphas == 0:\n",
    "        passes += 1\n",
    "    else:\n",
    "        passes = 0\n",
    "#     print(num_changed_alphas, passes, end=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.        ,  2.49930576,  1.00046283,  1.49884293])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alphas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3.9995371724799988"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# np.array([alphas[i] * ys[i] * Xs[i] for i in range(m)]).sum(axis=0)\n",
    "# most alphas should be 0\n",
    "ws = sum([alphas[i] * ys[i] * Xs[i] for i in range(m)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm_ws = ws / np.sqrt(ws.dot(ws))\n",
    "(norm_ws**2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "margin = np.min(np.abs(Xs.dot(ws) + b) / np.sqrt(ws.dot(ws)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 2.2)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGX6xvHvMym00CShLIKAdJAakLZIGoSEEIogLSBF\npCniroK4CCoqiv5EQVRQFFYUkV4TqmUtK6CIdAGRIkiTrrS8vz9mho0IZoAk7yTn+VwXF8mZM8mN\ncO68PnPmHDHGoJRSyjlctgMopZTKWlr8SinlMFr8SinlMFr8SinlMFr8SinlMFr8SinlMFr8Sinl\nMFr8SinlMFr8SinlMIG2A1xNSEiIqVy5su0YjrZnzx7Onj1LhQoVCAgIsB3H0bZt20bx4sUpWLCg\n7SjKj61bt+6IMSbMl339svgvXLjAfffdx/333287imMZYxg4cCDr168nOTmZAgUK2I7kWF9++SWJ\niYmMHz+e2NhY23GUnxKRn3ze1x+v1XPHHXeYkydPsnDhQmrUqGE7jmOlpqYycOBANmzYQHJyMvnz\n57cdybG++OILOnTowIYNGyhSpIjtOMoPicg6Y0y4L/v65Yw/V65cfP/999SoUYM9e/bYjuNYLpeL\n1157jRo1ahAbG8upU6dsR3KsRo0asWXLFooUKaLHhLppfln8AAUKFODcuXNERkYyefJk23Ecy1v+\nd9xxh5a/ZQUKFODgwYPUq1ePZcuW2Y6jsjG/LX5wr/yTk5N5+umntfwtcrlcTJw4kerVq9OyZUst\nf4uKFy/O3Llz6datm5a/umF+XfwA5cuXZ9WqVXz00UecO3fOdhzHcrlcvP7661SrVo24uDgtf4sa\nNWrE3Llzeeutt/DH1+iU//P74gd3+S9btowLFy6wYMEC23Ecy1v+VapU0fK3rHHjxsycOZMDBw7w\n8ccf246jsplsUfxehw8fZuDAgbz99tu2oziWy+XijTfeoHLlysTFxXH69GnbkRxtz549dOzYkRUr\nVtiOorKRbFX8ZcuWZdWqVYwaNYopU6bYjuNYLpeLN998U8vfDzRo0IA5c+bQpUsXLX/ls2xV/AAV\nKlRg5cqVXLhwwXYUR/OWf8WKFYmPj9fyt6hJkybMnj2bM2fO2I6isgm/fANXeHi4Wbt2bbr7LViw\ngKNHj9KzZ88sSKWuJjU1lb59+/LDDz+wePFiQkJCbEdytAkTJlClShWioqJsR1FZLNu/gctXlSpV\nYsSIEbz77ru2oziWy+Vi0qRJlC9fnvj4eF11WlajRg06d+7MqlWrbEdRfizbF//KlSv517/+xaxZ\ns2zHcSyXy8XkyZO5/fbbtfwta9q0KbNmzaJTp058/fXXtuMoP5WtRz1e27dvJ3/+/BQrVgyXK1v/\nLMvWUlNT6dOnD7t27WLx4sXky5fPdiTHWrt2LVWqVCFPnjx6TDiEY0Y9XhUrVqREiRJ06dKFqVOn\n2o7jWC6Xi7feeouyZcvSqlUrXflbFB4eTp48eWjatKme56/+JN3iF5FSIrJaRDaLyCYRGXyVfURE\nXhWRHSKyQUTqpHksVkS2eR4bltF/gLRGjRrF8OHDtfwt8pb/bbfdRkJCAmfPnrUdybFcLhejR4+m\nQ4cO7vI/dwrWTYXlI92/n9M34DlVuqMeESkBlDDGfCMi+YF1QBtjzOY0+8QBDwBxwJ3AK8aYO0Uk\nANgOxAD7gDVA57TPvZrrHfWktXXrVqKjo1mxYgV6Mxd7Ll26RO/evdmzZw+LFi0ib968tiM51scf\nf0y3zh3Z2CeAQrmBC2chKC+IC7rOgtsa2o6oMkCGjnqMMQeMMd94Pj4FbAFKXrFbIjDNuH0FFPL8\nwKgP7DDG7DLGnAdmePbNNJUrV2bjxo1UrlyZffv2Zea3Un8hICCAt99+m1KlStGqVStd+VvUrGFd\nd+kHnGXfUc/7LS6chfOnYfrdcE7fg+E01zXjF5EyQG3gv1c8VBLYm+bzfZ5t19qeqQoVKsRvv/1G\n48aNee+99zL726lrCAgIYMqUKdx666069rFp4xwK5Ya9J1Kp8+YZPv3p4v8eM6mwaY69bMoKn4tf\nREKA2cBDxpiTGR1ERPqKyFoRWXv48OGb/np58uRh6dKlPProo1r+FgUEBPDOO+9QsmRJWrdureVv\nw7FdcOEspQq6mHF3Hu6e+dv/yv/CWTi6y24+leV8Kn4RCcJd+tONMVdbHuwHSqX5/FbPtmtt/xNj\nzCRjTLgxJjwszKf7BaeratWqrFixgqlTp+olHizylv/f/vY3LX8bbinnnukDkWUD+aB9Hl5fe979\nWFBeKFLOYjhlgy9n9QjwNrDFGPN/19htAdDdc3ZPA+CEMeYA7hdzK4hIWREJBjp59s0yVatWZfny\n5fz+++8sXbo0K7+1SsNb/iVKlCAxMZHffvvNdiTnqN7O/UKuR1S5QD5on5d9J1P5Ys9FqNbOYjhl\ngy8r/sZAEhApIus9v+JEpJ+I9PPsswTYBewAJgMDAIwxF4FBQAruF4VnGmM2ZfQfwhcHDhygd+/e\nvP/++za+vcJd/u+++y7FihWjdevWWv5ZJVd+99k7wSGXV/4E5WXnqVy0mX2Rz77+1m4+leVyxDt3\nfbVp0yZiYmJ48cUX6dKlS4Z/feWbS5cu0b17dw4fPsz8+fPJkyeP7UjOcO60+4Xco7vc451q7Vjx\n2Vd06dKFOXPm0KRJE9sJ1U1w3Dt3fVWtWjWWLVum7yi1LCAggGnTphEWFqZjn6yUKwTqdIeYUe7f\nc4UQHR3N9OnTOX78uO10Kgs5asWf1uzZszl//jydO3fO1O+jru3ixYt0796do0ePMm/ePF35WzZu\n3Djq1atH48aNbUdRN0BX/D6oWLEiDz/8MDNmzLAdxbECAwOZNm0aRYoUoW3btvz++++2Izla1apV\nadu2LZ9//rntKCqTObb477jjDpYtW8aQIUNYuHCh7TiO5S3/woUL06ZNGy1/i5o3b857771H27Zt\nWb9+ve04KhM5dtTjtWnTJsLCwggNDdXL11p08eJFkpKS+PXXX5k3bx65c+e2HcmxvvrqK2rWrEmu\nXLn0mMhGdNRzHapVq0bRokXp0KEDM2fOtB3HsQIDA/n3v/9NoUKFdOxjWYMGDQgKCqJx48Z8+eWX\ntuOoTOD44vcaOXIkgwcP1vK3KDAwkPfee48CBQrQrl07LX+LAgMDGTlyJImJiVr+OZAWv0eNGjVI\nSUlh8ODB7Nql1y6xJTAwkOnTp5M/f34tf8tiY2OZNm0aHTt25PRpvYJnTqLFn0aNGjXYvHkz5cqV\nY//+q15SSGUBb/mHhITQvn17zp07ZzuSY8XGxrJhwwZCQkL0mMhBtPivULhwYU6fPs2dd96pN3C3\nyFv++fLlo127dlr+FhUuXJidO3dSu3Zt/vvfK6/IrrIjLf6rCAkJYfHixQwaNIjZs2fbjuNYQUFB\nTJ8+nbx582r5W3b77bfz7rvvkpCQoOWfA2jxX0PNmjVJTk7mjTfe4NKlS7bjOFZQUBDvv/8+efLk\n0bGPZXFxcbzzzjtMmDDBdhR1kxx/Hn96jDGcPHmSNWvWEB0dbTuOY124cIHOnTtz7tw5Zs2aRa5c\nuWxHcixjDHv27OHQoUPUq1fPdhzloefxZyARYe/evXTr1o25c+fajuNYQUFBfPDBBwQHB9OhQwdd\n+VskImzbto34+Hi+/vpr23HUDdDi90H16tVZunQp/fv31/K3KCgoiBkzZhAUFKTlb1nz5s2ZMmUK\nrVq1Ys2aNbbjqOukxe+j2rVrs2TJEo4dO2Y7iqN5yz8wMJAOHTpw/vx525Ecq1WrVrz99tscPXrU\ndhR1ndKd8YvIFKAVcMgYU/0qjz8CdPV8GghUAcKMMcdEZDdwCrgEXPR1/uRPM/6r+fDDD8mVKxdt\n2rSxHcWxzp8/T6dOnbh06RIfffQRwcHBtiM52ksvvUTTpk115m9RRs/43wVir/WgMWasMaaWMaYW\n8BjwiTEm7bI4wvO4T4Gyg/Lly3P//fczf/5821EcKzg4mBkzZuByuXTl7wfKly9Pq1at8OcFm/qf\ndIvfGPMp4Ot8ozPwwU0lygbq1q3LkiVL6Nu3L8uXL7cdx7GCg4P58MMPERE6duyo5W9RYmIikyZN\nIj4+ns2bN9uOo9KRYTN+EcmL+/8M0r7jyQArRGSdiPTNqO/lD+rWrUtKSgq1a9fGH0+JdYrg4ODL\nF9a75557tPwtSkxMZPbs2ZQrV06PCT+XkS/uJgCfXzHmaeIZAbUEBopI02s9WUT6ishaEVl7+PDh\nDIyVeWrVqkVoaCiJiYl6MxeLvOVvjNHyt6xJkya4XC4aNmzIN998YzuOuoaMLP5OXDHmMcbs9/x+\nCJgL1L/Wk40xk4wx4caY8LCwsAyMlfmeeOIJ+vTpo+Vvkbf8U1NT6dSpExcuXLAdybGCg4MZNmwY\nLVu21PL3UxlS/CJSELgLmJ9mWz4Rye/9GGgObMyI7+dvwsPDWbRoEffddx979+61HcexgoOD+eij\nj7h48SL33HOPlr9Fbdq04Y033qBdu3b89ttvtuOoK6Rb/CLyAfAlUElE9olIbxHpJyL90uzWFlhm\njDmTZlsx4D8i8h3wNbDYGJOckeH9Sb169di0aROlSpXiwIEDtuM4Vtry15W/XW3btmXdunXkyZNH\njwk/o9fqyWDHjx+nSpUqTJ48mVatWtmO41jnzp3j7rvvvnzaZ1BQkO1IjrVlyxaaNWtGcnIytWvX\nth0nx9Jr9VhUqFAh5s+fT69evVi8eLHtOI6VK1cuZs2axblz5+jcubOu/C2qUqUKEydOJDY2lm+/\n/dZ2HIUWf6aoX78+Cxcu5JVXXiE1NdV2HMfKlSsXs2fP5vfff6dLly5a/ha1b9+eiRMn8uqrr9qO\notBRT6YyxnDixAm+++477rrrLttxHOvcuXO0b9+ePHny8P777+vYxyJjDD/99BMnTpygZs2atuPk\nKDrq8RMiwq5du+jQoQNLly61HcexvCv/s2fP0rVrV135WyQifP/997Ro0YLvvvvOdhzH0uLPZHXq\n1GHBggX06NGD5OQce1KT3/OW/+nTp+natSsXL160HcmxEhISGD9+vJa/RVr8WaBBgwbMnz+fgwcP\n2o7iaLlz52bOnDla/n6gQ4cOjB8/nuzyLv2cRmf8WWz69OmEhobSokUL21Ec6/fff6dt27YUKFCA\n6dOnExgYaDuSo40dO5YWLVpQo0YN21GyNZ3x+7GyZcuSlJRESkqK7SiOlTt3bubOncvJkyfp1q2b\nrvwtK126NC1atOD777+3HcUxtPizWKNGjZg3bx5JSUl8+umntuM4lrf8jx8/TlJSkpa/Rffccw/j\nxo2jefPm/PDDD7bjOIKOeixZs2YNt99+O4ULF0ZEbMdxrN9//502bdpQuHBh/v3vf+vYx6KPP/6Y\nhg0bEhwcrMfEDdBRTzZQr149ChcuTHx8vN7MxaLcuXMzb948jh07Rvfu3XXlb1GzZs0wxtCwYUM2\nbsyR13P0G1r8FokIw4cPp2vXrlr+FnnL/8iRI/To0UPL36LcuXPz0EMP0bx5cy3/TKTFb1mTJk2Y\nM2cOSUlJerqnRXny5GH+/PkcPnyYHj16cOnSJduRHKtTp0689NJLJCYm6k11MonO+P3EkSNHCA0N\n5ZdffqFYsWK24zjWb7/9RmJiImFhYUybNo2AgADbkRzryJEjFClShEOHDukx4QOd8WdDoaGhHDly\nhJo1a7Jy5UrbcRzLu/I/dOiQrvwtCw0NZcOGDdSuXZtNmzbZjpOjaPH7kdDQUGbOnEmnTp20/C3y\nlv8vv/zCvffeq+VvUc2aNRk7diwxMTFs3rzZdpwcw5c7cE0RkUMictVXWkSkmYicEJH1nl9PpHks\nVkS2icgOERmWkcFzqqZNmzJ79mxefPFF/HEM5xR58+Zl/vz5HDhwQMvfsq5du/LCCy8wbtw421Fy\njHRn/CLSFDgNTDPGVL/K482AfxpjWl2xPQDYDsQA+4A1QGdjTLo/tp0447+SMYZff/2VrVu30qhR\nI9txHOvs2bMkJCRQsmRJ3nnnHZ35W2SMYdeuXZw/f54qVarYjuN3MnTGb4z5FDh2AznqAzuMMbuM\nMeeBGUDiDXwdRxIRtm/fTps2bVi9erXtOI6VN29eFi5cyP79++nVq5eu/C0SEdatW0d0dDRbtmyx\nHSdby6gZfyMR2SAiS0WkmmdbSWBvmn32ebYpHzVo0ICZM2dyzz338PHHH9uO41je8t+7dy+9e/fW\n8reoY8eOjBkzhujoaLZu3Wo7TraVEcX/DVDaGFMDGA/Mu5EvIiJ9RWStiKzVS7X+T7NmzZg5cyZ7\n9+5Nf2eVafLmzcuiRYvYs2ePlr9lSUlJjBkzRt/3chNuuviNMSeNMac9Hy8BgkQkFNgPlEqz662e\nbdf6OpOMMeHGmPCwsLCbjZWjNGvWjKSkJKZOnaorf4u8K/+ffvqJPn36aPlblJSURLNmzXj++ed1\n5X8Dbrr4RaS4eK6oJCL1PV/zKO4XcyuISFkRCQY6AQtu9vs5WalSpejQoQOffPKJ7SiOlS9fPhYt\nWsTu3bu57777SE1NtR3J0YoVK0Z0dDTbtm2zHSVb8eV0zg+AL4FKIrJPRHqLSD8R6efZ5W5go4h8\nB7wKdDJuF4FBQAqwBZhpjNF3YdyEyMhIPvzwQzp06MB///tf23Ecy1v+P/74I3369NHyt+jee+9l\n9OjRREVF8dNPP9mOk23oJRuyoS+++IIqVapQqFAhvXytRWfOnCE+Pp5y5crx1ltv4XLp+yFtWbFi\nBU2bNiUoKMixx4ResiGHa9SoEYUKFSI2NlZv5mJRvnz5WLx4MTt37tSxj2XR0dFcuHCBRo0asX37\ndttx/J4WfzYlIjz66KO0b9+ezz77zHYcx/KW/44dO+jbt6+Wv0X58uXjvvvuIyoqSss/HVr82VhU\nVBQffPABHTp04OjRo7bjOFZISAiLFy9m+/btWv6W9erViyeffJLWrVvrfRX+gs74c4BDhw5RtGhR\nDh8+jJ4Ka8/p06eJi4ujUqVKvPnmmzrzt+jQoUOEhYVx5MgRxxwTOuN3mKJFi3LgwAGqV6/Of/7z\nH9txHCskJIQlS5awdetW+vXrpyt/i4oWLcqaNWuoU6eO3sD9KrT4c4gSJUrw73//m7Zt22r5WxQS\nEsLSpUvZsmWLlr9l9evXZ+TIkURGRrJjxw7bcfyKFn8O0rx5c6ZPn86YMWP0ks4WeVf+mzdvpn//\n/lr+FvXp04cnnniCl19+2XYUv6Iz/hzIGMPRo0f58ccfqVevnu04jnXq1ClatmxJ9erVmThxos78\nLTLG8MMPPxAQEMDtt99uO06m0Bm/w4kImzZtIj4+ni+++MJ2HMfKnz8/S5cu5fvvv2fgwIG68rdI\nRPjqq6+IjIxk586dtuNYp8WfQ911111MmzaNNm3aaPlblD9/fpKTk9mwYYOWv2Xdu3dn+PDhWv5o\n8edosbGxTJs2Ta9hYpl35f/dd98xaNAgff3Fovvvv58RI0Zw4MAB21Gs0hm/Q0yZMoUqVarQsGFD\n21Ec6+TJk8TGxlKrVi1ee+01x15Txl+MGTOGjh07Uq5cOdtRMoTO+NWflChRgsTERL766ivbURyr\nQIECJCcns379egYOHKgrf8sKFChAZGQku3btsh0ly2nxO0TLli2ZOnUqrVu35ptvvrEdx7G85f/t\nt9/q2MeyAQMG8OijjxIZGcn+/de8R1SOpKMeh/nkk0+oVasWBQoU0FGDRSdOnCA2Npa6desyfvx4\n/buwaOnSpURHRxMYGJit/x501KOu6a677iJ//vw0b95cb+ZiUcGCBUlOTmbt2rU88MADuvK3qGXL\nlvz22280btyY3bt3246TJXy5A9cUETkkIhuv8XhXEdkgIt+LyBciUjPNY7s929eLiC7h/YTL5eKh\nhx4iISFBy9+iggULkpKSwtq1a3nwwQe1/C0qUKAAXbt2pVmzZo4of19W/O8CsX/x+I/AXcaYO4Cn\ngUlXPB5hjKnl6/+CqKwRHx/PlClTaNu2LSdOnLAdx7G85f/1118zePBgLX+LBg4cyD//+U9at27N\npUuXbMfJVD7N+EWkDLDIGFM9nf0KAxuNMSU9n+8Gwo0xR64nlM74s87BgwcpXrw4R44cITQ01HYc\nxzp+/DgtWrSgQYMGjBs3LlvPmrO7gwcPUqxYMY4dO0aRIkVsx/GZzRl/b2Bpms8NsEJE1olI3796\nooj0FZG1IrL28OHDGRxLXUvx4sXZu3cv1atXZ82aNbbjOFahQoVISUnhyy+/5KGHHtKVv0XFixfn\nP//5D+Hh4Tn2zY8ZVvwiEoG7+Iem2dzEGFMLaAkMFJGm13q+MWaSMSbcGBPulBsn+ItSpUoxefJk\n4uPjtfwtKlSoEMuWLePLL79kyJAhWv4W/f3vf2fIkCFERETkyPLPkOIXkRrAW0CiMebyPQCNMfs9\nvx8C5gL1M+L7qYyXkJDAW2+9xTPPPGM7iqN5y//zzz/n4Ycf1vK36MEHH2Tw4MGMGzfOdpQMd9Mz\nfhEpDawCuhtjvkizPR/gMsac8ny8HHjKGJOc3vfTGb89xhgOHz7Mzz//TK1atWzHcazjx48TExND\nkyZN+L//+z+d+VtkjGHbtm3kzZuX0qVL245zTRk64xeRD4AvgUoisk9EeotIPxHp59nlCaAIMPGK\n0zaLAf8Rke+Ar4HFvpS+sktEWL9+PS1atGDdunW24ziWd+X/2Wef8Y9//ENX/haJCJ988gkRERHs\n2bPHdpwMoe/cVVc1f/58+vbty9KlS6lTp47tOI7166+/EhMTw1133cWLL76oK3+LXn75ZSZMmMDH\nH39MqVKlbMf5k+tZ8QdmdhiVPSUmJgKwc+dOLX6LChcuzPLly4mOjuaRRx5h7NixWv6WDBkyhKCg\nIPbv3++XxX89dMWv0jV58mTq1q2rPwAs+vXXX4mOjiYiIkLL3zJjDM899xxJSUl+9QNAr9WjMlRo\naChxcXF8++23tqM4lnflv2rVKh599FGd+VskIuTKlYuIiAj27t1rO84N0eJX6Wrbti2vv/46LVu2\nZOPGq16ySWWBW265hRUrVrBy5UqGDh2q5W/RP/7xD/r3709ERAS//PKL7TjXTWf8yidt27alQIEC\nlC5dGmOMjhos8ZZ/dHQ0Q4cO5fnnn9e/C0v+8Y9/ULFiRYoUKZLtjgld8SufRUVFkTdvXqKjo1m/\nfr3tOI7lLf/ly5czbNgwXflblJCQwMmTJ2natCn79u2zHcdnWvzqugQGBjJgwABiY2O1/C3ylv+y\nZct47LHHtPwtuuWWW0hISCAiIiLb3MlLRz3qurVv3x5wr3a2bNlCSEiI5UTOVKRIEVasWEFUVBQi\nwrPPPputxg05yaOPPgq4T4Nes2aN3/896Omc6oYdOHCAEiVKcOzYMW655RbbcRzr6NGjREVF0bJl\nSy1/yw4cOEDx4sU5fvw4hQsXztLvradzqixRokQJdu7cSbVq1diwYYPtOI7lXfkvWbKExx9/XMc+\nFpUoUYKVK1fSoEEDfv75Z9txrkmLX92U22+/nVdeeYXmzZtr+VsUGhrKypUrWbx4Mf/617+0/C2K\njo6mZ8+eNGvWzG/LX2f86qZ17NgRYwxPP/00H330ke04juUt/6ioKABGjx6tYx9LvGdbvfLKKzz/\n/PO24/yJzvhVhklNTeXQoUMcPXqUatWq2Y7jWEeOHCEyMpLWrVvz9NNPa/lblJqaytatWylcuDAl\nSpTI1O+lM35lhcvlYs2aNURHR+s7fC3yrvwXLFjAE088oWMfi1wuF8uXLyciIoIDBw7YjnOZFr/K\nUAkJCbz88ss0b95cy9+isLAwVq5cybx58xg5cqSWv0WDBw8mKSmJyMhIvyn/dGf8IjIFaAUcusYd\nuAR4BYgDzgL3GmO+8TwW63ksAHjLGDMmA7MrP9WpUycAtm/fzj5CaVGt2B/GDcYYUjb98qftKmOF\nhYWxatUqIiMjAXjyySf1v7cljz/+OMHBwezfv9898jl3CjbOgWO74JZyUL0d5MqfZXnSnfF7bpB+\nGph2jeKPAx7AXfx3Aq8YY+4UkQBgOxAD7APWAJ2NMZvTC6Uz/pwheeNBuv5zNJ0SYpgwoDUigjGG\npxZt5p3Pd/NGt7rEVi9uO2aOd+jQISIjI2nfvj1PPvmk7TiOZozhuaH96RUwj+IhAhfOQlBeEBd0\nnQW3Nbzhr52hM35jzKfAsb/YJRH3DwVjjPkKKCQiJXDfWH2HMWaXMeY8MMOzr3KIFtWK0ax6ad4c\n1otBExf8ofR7Ni5Di2rFbEd0hKJFi7Jq1Spmz57NqFGjbMdxNDl/mkvfvE/kW4c5+Otp98YLZ+H8\naZh+N5w7nSU5MuJ0zpJA2otS7/Nsu9r2OzPg+6lsQkSYNfYR7gbeHNaLed89R9AtJenZuAxPtKqq\nY4cs5C1/79hHfwBYsnEOIyLywaXzRE49y2c981Ikr2f9bVJh0xyo0z3TY/jNefwi0hfoC/j1nezV\n9fGWf/E9vxMQcgvGGC19S7zlHxERgYgwcuRI25Gc59guuHCWEXflonpRF4Vyy/8u6XzhLBzdlSUx\nMuKsnv1A2vuP3erZdq3tV2WMmWSMCTfGhIeFhWVALOUPvOOdPOXqIgGBHPrwcQZNXKhnmVjiLf+Z\nM2fy1FNP2Y7jPLeUc8/0gbZVgjj6myFy2ll+OZ3q3l6kXJbEyIjiXwB0F7cGwAljzAHcL+ZWEJGy\nIhIMdPLsqxziypn+7hcSadmuM28O68kDr2v521KsWDFWrVrFjBkztPyzWvV27hdyPYrmc3HXbQHu\n8j9joFq7LInhy+mcHwDNgFAR2QeMBIIAjDFvAEtwn9GzA/fpnD09j10UkUFACu7TOacYYzZlwp9B\n+amUTb9cLn3veGf2S0NpD7zx2H3E1P+KxPCytmM6UrFixVi9evXlsc+IESNsR3KGXPndZ+9Mv9s9\n079wllExt0DAadotK8p/nslHVgxB9ZINKtNc63x9YwzTV62na2QtTpw4QaFChSymdLaDBw8SERFB\n165d+ddJnKo6AAAXJ0lEQVS//mU7jnOcO+1+IffoLvd4p1o79h85QYkSJTh16hQFCxa87i+pl2xQ\nfkFEiK1e/E8v5IoI3aJqs3XrVqpXr87WrVstJVTFixdn9erVTJ8+ndGjR9uO4xy5Qtxn78SMcv+e\nK4SSJUuyZMkSmjRpwqFDhzL122vxK2uqVKnCM888Q1RUlJa/RWnL/5lnnrEdx9Hi4+Np06YNUVFR\nHD58ONO+j9+czqmcqUePHu4XgZ96ivfff992HMcqXrz4H071HD58uO1IjiQil19wHz9+fKa9+K4z\nfuUXUlNTOXjwIKdPn6ZixYq24zjWgQMHiIiIoEePHjz22GO24ziWMQZjDFu2bKFo0aL4coq7zvhV\ntuNyufj888+JjIxk27ZttuM4VokSJVi9ejVTp07lueeesx3HsUQEl8vF/PnziYqK4siRIxn69XXU\no/xGhw4dOHXqFFFRUaxatUpX/pZ4y79Zs2aICMOGDbMdybEee+wxzpw5Q1RUFCtXriQ0NDRDvq4W\nv/IrvXr1QkTYvn27Fr9F3vL3zvyHDh1qO5IjiQijR48mT5487N27N8OKX2f8ym+9/vrrREdHU6FC\nBdtRHOvnn38mIiKCXr16aflbZoxhzJgx9O3blyJFivzpcZ3xqxwhODiYyMhIfvjhB9tRHOtvf/sb\nq1ev5u233+aFF16wHcfxTp48SVRUFEePHr2pr6OjHuW3evfuDUBkZCSffvopZcvq5R1s8Ja/d+zz\nyCOP2I7kSCLCs88+izGGqKgoPvnkkxt6hy9o8Ss/17t3b4oWLUrRokVtR3G0kiVLXi5/QMvfEhHh\nueeeo27duoSEhNzw19FRj/J7CQkJBAQEEBUVxY4dO2zHcSxv+U+aNIkXX3zRdhzHEhE6dOjAwYMH\niYmJuaGxjxa/yhZy587NPffcQ2RkpJa/Rd7yf/PNN3nppZdsx3G0v/3tb9SuXZuYmBiOHfuru+P+\nmY56VLbRt29fAGJjY9m8eTPBwcGWEznTrbfe+ofz/B9++GHbkRxJRHj++ed59NFH6dix4/U9V0/n\nVNnN3r17KVWqFCdOnLjhF7fUzdu7dy8REREMGDBAy98iYwz79u2jdOnSejqnyrlKlSrF+vXrqVGj\nBjt37rQdx7FKlSrF6tWrmThxIi+//LLtOI4lIpQqVSr9HdPwqfhFJFZEtonIDhH50/u3ReQREVnv\n+bVRRC6JyC2ex3aLyPeex3QZrzJErVq1GDZsGJGRkezalTU3qFZ/5i3/CRMmaPlnI77cejEAeA2I\nAfYBa0RkgTFms3cfY8xYYKxn/wRgiDEm7asNEcaYjL3KkHK8/v37Y4zhySefZOrUqbbjOJa3/L3n\n+T/00EO2I6l0+PLibn1ghzFmF4CIzAASgc3X2L8z8EHGxFPqrw0YMID777+fffv2ceHCBX2TlyWl\nS5f+w3n+Wv7+zZdRT0lgb5rP93m2/YmI5AVigdlpNhtghYisE5G+NxpUqWsJCAi4XDo//vij7TiO\n5S3/V199lVdeecV2HPUXMvp0zgTg8yvGPE2MMftFpCiwXES2GmM+vfKJnh8KfcH9D0ip65GUlMTJ\nkyeJiIhg9erVuvK3JO3KX0R48MEHbUdSV+FL8e8H0r5kfKtn29V04ooxjzFmv+f3QyIyF/fo6E/F\nb4yZBEwC9+mcPuRS6g8GDhx4+ZLOWvz23HbbbX8o/wceeMB2JHUFX4p/DVBBRMriLvxOQJcrdxKR\ngsBdQLc02/IBLmPMKc/HzYHMuYmkUrhn/gATJ04kLi6OMmXK2A3kUGnLH9Dy9zPpFr8x5qKIDAJS\ngABgijFmk4j08zz+hmfXtsAyY8yZNE8vBswVEe/3et8Yk5yRfwClriY1NfXy2EfL347bbrvtDzdw\nHzRokO1IysOnGb8xZgmw5Iptb1zx+bvAu1ds2wXUvKmESt2AQYMGYYwhIiKCzz77jFtvvdV2JEcq\nU6bMH8Y+AwcOtB1JodfqUTnYAw88QMmSJa96tyKVddKWP6Dl7wf0kg0qR2vXrh2pqanExMTw008/\n2Y7jWN7yHzt2LBMnTrQdx/G0+FWOly9fPlq1akVERAR79uyxHcexvOX/wgsv8Prrr9uO42g66lGO\nMHjwYMB9SecNGzYQGKj/9G0oW7Ysq1atIjIyEnBfdkNlPf3Xrxxj8ODBtGnThsDAQE6ePEmBAgVs\nR3KkcuXK/eFsn379+tmO5Dg66lGOctttt/H1119Tp04d9u7dm/4TVKYoV64cq1ev5rnnnuPNN9+0\nHcdxtPiV49SvX58BAwYQERGh5W+Rd+X/7LPPavlnMR31KEfy3jHq6aefZtKkSZbTONftt99+eeYv\nIpdvr6kylxa/cqyHH36YS5cusWfPHlwul77JyxJv+XvP89fyz3w66lGOFhAQQHJyMhEREezbt892\nHMfylv/o0aOZPHmy7Tg5nq74leP17duXEydOXL62j6787Shfvvwfxj59+vSxHSnH0uJXCnjkkUcI\nCAhg27ZtWvwWpS1/QMs/k2jxK+XhfcH3tddeo02bNpQsedUbzalMVr58eVauXHl55d+7d2/bkXIc\nnfErdYXTp08TERHB/v3Xut+QymwVKlRg1apVjBo1iilTptiOk+Poil+pKwwdOhSAyMhIPvnkE4oX\nL245kTN5y9879unVq5flRDlHjit+Ywwpm36hRbVieG4A85fblbqaoUOHUqZMGQoVKmQ7iqNVqFCB\nlStXEhUVhYjQs2dP25FyBJ9GPSISKyLbRGSHiAy7yuPNROSEiKz3/HrC1+dmtJRNv9DvvXU8tWgz\nxrhv3WuM4alFm+n33jpSNv2S2RFUDnHPPfdw/vx5YmNj+fnnn23HcayKFSuycuVKRowYwbvvvms7\nTo6Q7opfRAKA14AYYB+wRkQWGGM2X7HrZ8aYVjf43AzToloxejYuwzuf7wbgiVZVeWrRZt75fDc9\nG5ehRbVimfWtVQ5UoEABmjZtSmRkJKtXr6ZEiRK2IzmSt/yjoqIAuPfee+0GyuZ8GfXUB3Z4bqOI\niMwAEgFfyvtmnntDRIQnWlUF4J3Pd1/+AdCzcRmeaFVVxzzqug0fPhyAuLg41q5dS0BAgOVEzlSp\nUqU/jH169OhhO1K25cuopySQ9kpW+zzbrtRIRDaIyFIRqXadz0VE+orIWhFZe/jwYR9iXVva8vfS\n0lc3Y/jw4cyePZuAgABOnz5tO45jecv/8ccfZ+rUqbbjZFsZdTrnN0BpY0wNYDww73q/gDFmkjEm\n3BgTHhYWdlNhvDP9tNLO/JW6EeXKlePTTz+lXr16HDx40HYcx6pUqRIrVqzg8ccfZ9q0abbjZEu+\nFP9+oFSaz2/1bLvMGHPSGHPa8/ESIEhEQn15bkbzlr53pv/jc3GXZ/5a/upmNW3alC5duhAREaHl\nb1HlypVZsWIFjz32mJb/DfBlxr8GqCAiZXGXdiegS9odRKQ48IsxxohIfdw/UI4Cx9N7bkZL2fTL\n5dL3jnfSzvzvLFuE2Op6Xra6cSNGjADgmWeeYfz48ZbTOFflypX/MPNPSkqyHSnbSLf4jTEXRWQQ\nkAIEAFOMMZtEpJ/n8TeAu4H+InIR+A3oZNxL66s+N5P+LID7rJ43utX9w/n63vK/s2wRPatHZYgR\nI0Zw8eJFdu/eTZ48eShWTP9d2eBd+UdHRyMidOvWzXakbEH8cfQRHh5u1q5dazuGUul65ZVXmDRp\nEqtWrdLyt2jz5s3ExMTw/PPPO7b8RWSdMSbcl31z3Dt3lcpKgwcP5vjx40RGRmr5W1S1alWWL19+\neeXftWtX25H8mha/Ujdp5MiRBAYGsnXrVi1+i6pWrfqHsU+XLpn6cmK2psWvVAZ4/PHHMcYwYcIE\nOnbsSNGiRW1HciTvyj8mJgZAy/8a9LLMSmWgw4cPExkZyaFDh2xHcaxq1aqxfPly/vnPf/L+++/b\njuOXdMWvVAYREUaNGoUxhqioKFavXk1oaKjtWI7kLf+YmBhEhM6dO9uO5Fe0+JXKQCLCk08+SaVK\nlcifP7/tOI5WrVo1li1bRvPmzREROnXqZDuS39BRj1IZzHtWyZkzZ4iPj+dmrz2lblz16tVZtmwZ\nQ4YMYcaMGbbj+A0tfqUySeHChalduzZRUVFa/hZVr16d5cuXM2TIED788EPbcfyCjnqUyiQiwtNP\nP40xhoSEBL788ku9Qqwl3pW/d+zTsWNH25Gs0uJXKhOJCKNHj758y8AzZ86QL18+y6mc6Y477iAl\nJYUWLVoAOLr8ddSjVCYTEcqXL8/y5ctp1KgRR44csR3JsWrUqEFKSgqDBw/mo48+sh3HGi1+pbJI\nTEwMcXFxREdHc/ToUdtxHMtb/g8++KBjy1+LX6ksIiI8++yztGzZkjFjxtiO42g1atQgOTmZBx98\nkFmzZtmOk+V0xq9UFvKW/6VLl9i1axcFCxakSJEitmM5Us2aNUlOTr4887/77rstJ8o6uuJXKouJ\nCIGBgcycOVPHPpbVrFmTlJQUBg0axOzZs23HyTI+Fb+IxIrINhHZISLDrvJ4V8+N1r8XkS9EpGaa\nx3Z7tq8XEb3IvlIeQ4cOJSYmhpiYGI4dO2Y7jmN5V/4DBw5kzpw5tuNkiXSLX0QCgNeAlkBVoLOI\nVL1itx+Bu4wxdwBPA5OueDzCGFPL15sEKOUEIsLzzz9PXFwcW7dutR3H0WrVqsXSpUsZMGCAI8rf\nlxV/fWCHMWaXMeY8MANITLuDMeYLY8yvnk+/wn1TdaVUOrzn+Tds2JAJEyboyt+i2rVrXy7/uXPn\n2o6TqXwp/pLA3jSf7/Nsu5bewNI0nxtghYisE5G+1x9RKWfYvXs30dHRWv4Wecu/f//+Obr8M/TF\nXRGJwF38Q9NsbmKMqYV7VDRQRJpe47l9RWStiKzV65oopxERxo4dS0REBDExMRw/ftx2JMeqXbs2\nS5YsoX///sybN892nEzhS/HvB0ql+fxWz7Y/EJEawFtAojHm8mkKxpj9nt8PAXNxj47+xBgzyRgT\nbowJDwsL8/1PoFQOISK8+OKLPPDAA+TNm9d2HEerU6cOS5Ys4f7778+R5e9L8a8BKohIWREJBjoB\nC9LuICKlgTlAkjFme5rt+UQkv/djoDmwMaPCK5XTiAj33nsvx48fJzExkV9//TX9J6lMUadOHZYu\nXcr999/P/PnzbcfJUOm+gcsYc1FEBgEpQAAwxRizSUT6eR5/A3gCKAJM9Fx98KLnDJ5iwFzPtkDg\nfWNMcqb8SZTKQcLCwihXrhwxMTEsX76cwoUL247kSN6Vf1xcHCJC69atbUfKEGKMsZ3hT8LDw83a\ntXrKv3I2YwxDhgxh/fr1rF69Wi/pbNHatWuJj49n8uTJflv+IrLO11Pm9ZINSvkpEeHll1/mhx9+\nAODs2bM6+7ckPDycxYsXEx8fj4iQkJBgO9JN0Us2KOXHRISKFSuyaNEimjVrpmf7WOQt/z59+rBw\n4ULbcW6KFr9S2UCrVq1o2LAhzZs31/K3KDw8nEWLFtGnTx8WLVpkO84N0+JXKhsQEcaNG0eDBg0Y\nO3as7TiOVq9ePRYtWkTv3r2zbfnri7tKZSPGGC5dusRPP/1EaGgoBQsWtB3Jsb7++msSEhKYMmUK\n8fHxtuNc14u7uuJXKhvxXtJ52rRptGjRghMnTtiO5Fj169dn4cKF9OzZkyVLltiOc120+JXKhkaN\nGkW9evW0/C3zlv+9996brcpfi1+pbEhEePXVV/n73//Otm3bbMdxtDvvvDPblb8Wv1LZlPfCbvXq\n1eO1117j5MmTtiM51p133smCBQu49957Wbp0afpPsEyLX6kcYPPmzbRo0ULL36IGDRqwYMECevTo\nQXKyf1+ZRotfqWxORJgwYQK1a9cmNjaWU6dO2Y7kWN7y7969u1+Xvxa/UjmAt/x79epF7ty5bcdx\ntAYNGjB//ny6d+9OSkqK7ThXpcWvVA7hcrno06cPhw8f5u6779aVv0UNGzZk3rx5JCUl+WX5a/Er\nlcMUL16c0NBQHftY1qhRo8vlv2zZMttx/kCLX6kcxuVyMXHiRKpXr07Hjh1tx3E0b/l369bNr8pf\ni1+pHMjlcvH6668zbtw4jDH8/vvvtiM5VqNGjZg7dy7dunVj+fLltuMAPha/iMSKyDYR2SEiw67y\nuIjIq57HN4hIHV+fq5TKHC6Xi0qVKjFr1ixiYmI4ffq07UiO1bhxY+bMmUPXrl39ovzTLX4RCQBe\nA1oCVYHOIlL1it1aAhU8v/oCr1/Hc5VSmah9+/ZUrlyZli1bavlb1KRJk8vlv2LFCqtZfFnx1wd2\nGGN2GWPOAzOAxCv2SQSmGbevgEIiUsLH5yqlMpHL5eLNN9+kcuXKvPzyy7bjOJq3/Lt06WK1/H0p\n/pLA3jSf7/Ns82UfX56rlMpk3vJ/7LHH2LFjh678LWrSpAmzZ8+mS5curFy50koGv7nnroj0xT0m\nAjgnIhtt5vkLocAR2yGuwZ+zgX/n8+ds4N/5NNsNio6Ozsh8t/m6oy/Fvx8olebzWz3bfNknyIfn\nAmCMmQRMAhCRtb7eUCCrabYb58/5/Dkb+Hc+zXbjbOXzZdSzBqggImVFJBjoBCy4Yp8FQHfP2T0N\ngBPGmAM+PlcppVQWSnfFb4y5KCKDgBQgAJhijNkkIv08j78BLAHigB3AWaDnXz03U/4kSimlfOLT\njN8YswR3uafd9kaajw0w0Nfn+mDSde6flTTbjfPnfP6cDfw7n2a7cVby+eXN1pVSSmUevWSDUko5\njLXi9/VSDiJST0Quisjd/pZPRJqJyHoR2SQin/hLNhEpKCILReQ7T7aeWZhtiogcutbpuH91eQ8/\nyNbVk+l7EflCRGpmVTZf8qXZL8uPCV+yWTwe0vt7tXk8lBKR1SKy2fO9B19ln6w/JowxWf4L9wu9\nO4FyQDDwHVD1Gvutwv0awd3+lA8oBGwGSns+L+pH2YYDz3s+DgOOAcFZlK8pUAfYeI3H44ClgAAN\ngP9m4d9retkaAYU9H7fMymy+5Evz92/jmEjvv52V48HHbDaPhxJAHc/H+YHtVzles/yYsLXi9/VS\nDg8As4FDWRkO3/J1AeYYY/YAGGOyKqMv2QyQX0QECMH9D/1iVoQzxnzq+X7Xcq3Le1jPZoz5whjz\nq+fTr3C/7yTL+PDfDiwdEz5ks3U8+JLN5vFwwBjzjefjU8AW/nz1giw/JmwVf7qXchCRkkBbPBd8\ny2K+XGqiIlBYRD4WkXUi0t2Psk0AqgA/A98Dg40xqVkTL13Z5TIevXGvwvyG5WMiPbaOB1/4xfEg\nImWA2sB/r3goy48Jv7lkw1WMA4YaY1LdP6j9TiBQF4gC8gBfishXxpjtdmMB0AJYD0QCtwPLReQz\nY8xJu7GyBxGJwF38TWxnuYI/HxN6PPwFEQnB/X9qD/nDcWir+H25DEQ4MMPzDzwUiBORi8aYeX6S\nbx9w1BhzBjgjIp8CNXHP8Gxn6wmMMe4B4g4R+RGoDHydydl84Ut+a0SkBvAW0NIYc9R2nivYPCbS\nY+t48IXV40FEgnCX/nRjzJyr7JLlx4StUU+6l3IwxpQ1xpQxxpQBZgEDsvAfuC+XmpgPNBGRQBHJ\nC9yJe37nD9n24F55ISLFgErArizI5otrXd7DOhEpDcwBkvxkpfoHlo+J9Ng6Hnxh7XjwvK7wNrDF\nGPN/19gty48JKyt+49tlIKzxJZ8xZouIJAMbgFTgLWNMpl9R1Mf/dk8D74rI97jPFBhqjMmSKxSK\nyAdAMyBURPYBI3FfrM+b7aqX9/CTbE8ARYCJnlX1RZOFF9DyIZ816WWzdTz4kg2LxwPQGEgCvheR\n9Z5tw4HSafJl+TGh79xVSimH0XfuKqWUw2jxK6WUw2jxK6WUw2jxK6WUw2jxK6WUw2jxK6WUw2jx\nK6WUw2jxK6WUw/w/OeZ/k4/PUsUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ff09518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_x1s = np.array([0, 3])\n",
    "_x2s = (-b - ws[0] * _x1s) / ws[1]\n",
    "_xs = np.array([_x1s, _x2s]).T\n",
    "assert (_xs[1] - _xs[0]).dot(ws) < 1e-12\n",
    "\n",
    "# decision boundary\n",
    "plt.plot(*_xs.T, lw=1, color='black')\n",
    "# margin\n",
    "plt.plot(*(_xs - norm_ws * margin).T, color='black', lw=1, linestyle='--')\n",
    "plt.plot(*(_xs + norm_ws * margin).T, color='black', lw=1, linestyle='--')\n",
    "\n",
    "plt.scatter(*Xs[:2].T, marker='x', s=50)\n",
    "plt.scatter(*Xs[2:].T, marker='o', s=50)\n",
    "\n",
    "plt.xlim(0.4, 2.1)\n",
    "plt.ylim(0, 2.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When I first see this margin, it looks so beautiful, and feels so good it works! 🍺"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
